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ABSTRACT 



cn 



Aims. We formulate an improved time series analysis method for the analysis of photometry of active stars. This new Continuous 

Period Search (CPS) method is applied to 12 years of V band photometry of the young solar analogue HD 1 16956 (NQ UMa). 

Ql ' Methods. The new method is developed from the previous Three Stage Period Analysis (TSPA) method. Our improvements are 

^^ I the use of a sliding window in choosing the modelled datasets, a criterion applied to select the best model for each dataset and the 

• , computation of the time scale of change of the light curve. We test the performance of CPS with simulated and real data. 

i-C ' Results. The CPS has a much improved time resolution which allows us to better investigate fast evolution of stellar light curves. We 

t^' can also separate between the cases when the data is best described by periodic (i.e. rotational modulation of brightness) and aperiodic 

(e.g. constant brightness) models. We find, however, that the performance of the CPS has certain limitations. It does not determine 

i^ , the correct model complexity in all cases, especially when the underlying light curve is constant and the number of observations too 

"*— * ' small. Also the sensitivity in detecting two close light curve minima is limited and it has a certain amount of intrinsic instability in 

_j ' its period estimation. Using the CPS, we find persistent active longitudes in the star HD 1 16956 and a "flip-flop" event that occurred 

during the year 1999. Assuming that the surface differential rotation of the star causes observable period variations in the stellar light 

curve, we determine the differential rotation coefficient to be \k\ > 0.1 1. The mean timescale of change of the light curve during the 

Cn ' whole 12 year observing period was Tc = 44. 1 d, which is of the same order as the predicted convective turnover time of the star. We 

^ ' also investigate the presence of activity cycles on the star, but do not find any conclusive evidence supporting them. 

O ! 

r^v Key words. Methods: data analysis. Stars: activity, starspots, individual: HD 116956 

o' 
^, 

r^^ 1. Introduction The second improvement is to test several models of different 

(^ order K, instead of just one, as possible descriptions of the data. 

(3 One of the most common methods to search for periodicity in The best model can then be selected with a suitable criterion. By 

T-H . an astronomical time series is the Lomb-Scargle periodogram including a constant brightness model into the set of tested pos- 

L; I (Lomb 1976; Scargle 1982). The method is, however, not suit- sible models, we can investigate whether a periodic function is 

. ^ ' able for analysing data where the periodic variations do not fit an appropriate description of the data at all. This is of particular 

^ to a simple single harmonic model. A more general period anal- interest when analysing low ampHtude stellar fight curves. 

5^ ysis metfiod, tfie Tfiree Stage Period Analysis (TSPA), was for- As tfie tfiird improvement to tfie TSPA, we define the pa- 

mulated by Jetsu & Pelt (1999, hereafter Paper 1). This method rameter Tq as the approximate time scale in which the shape of 

utifises a higher order harmonic model: a Kth order Fourier se- the analysed light curve undergoes significant changes. It char- 

ries. It can thus model more accurately periodic data with a com- acterises the time interval after which any particular light curve 

plex form. The main motivation in developing the TSPA was to model does not adequately describe the subsequent fight curves 

formulate a suitable method for analysing photometry of active any more. This parameter is, however, only an approximation, 

stars. as it depends also on the quality and quantity of the photometiic 

The TSPA method can still be improved. In this paper we data. Nevertheless, when analysing photometry of active stars, 

discuss and implement three such improvements. First of all, the parameter may tell us how fast the spot distribution on the 

TSPA is not flexible enough for studying fast evolution of the visible stellar surface is changing. 

fight curve, as it divides the photometry into separate datasets We formulate the new Continuous Period Search (CPS) 

that do not overlap. A better time resolution can be achieved by method in Sect. 2 by incorporating the above three improve- 

choosing the datasets with a sliding window and allowing tfie ments to tfie TSPA, and use simulated data to test tfie perfor- 

adjacent datasets to overlap. Tfiis produces a sequence of analy- mance of tfie CPS in Sect. 3. Finally, we apply tfie CPS to real 

sis results akin to video view in contrast to tfie set of snapsfiots pfiotometric observations of tfie spot activity on tfie surface of 

provided by tfie TSPA. tfie nearby young solar analogue HD 1 16956 (NQ UMa) in Sect. 

4. 
HD 1 16956 fias a spectral type of G9V, an effective temper- 



ed 



* The analysed photometry and numerical results of the analysis are ature of T^g = 5 170 K, a projected rotation velocity of v sin / = 
both published electronically at the CDS. 5.6 km s"' and a fitfiium abundance of log A^(Li) - 1.42 + 0.12 



J. Lehtinen et al.: The continuous period search method and its application to the young solar analogue HD 1 16956 



(Gaidos 1998). From its space motion, Gaidos et al. (2000) iden- 
tified the star as a Local Association member. This would sug- 
gest an age between 20 and 150 Myr (Montes et al. 2001). On 
the other hand, Gaidos et al. (2000) noted that none of the Local 
Association stars in their sample display characteristics of such 
a young age and rather suggest that they have similar ages to the 
Ursa Major group stars (300 Myr). Nevertheless, HD 116956 
must be considerably younger than the Sun, as also indicated by 
its lithium content. The radial velocity of HD 1 16956 is constant 
-12.3 km s"', i.e. it seems to be a single star (Gaidos et al. 2000). 
HD 116956 displays numerous signs of magnetic activity, 
such as rotational modulation of its brightness caused by large 
spots on its surface. The amplitude of this variability is typically 
0.05 mag. Gaidos et al. (2000) reported a preliminary rotation 
period estimate of P - 7.80 ± 0.02 days. Another activity in- 
dicator is the logarithmic ratio of X-ray and bolometric lumi- 
nosities Rx = -4.48, which indicates that the star has an ac- 
tive corona (Gaidos 1998). Finally, the logarithmic ratio of Ca II 
emission and bolometric luminosity, R'^^ = -4.447, means that 
HD 1 16956 is chromospherically active (Gray et al. 2003). 

2. The CPS method 

The light curves of active stars can undergo rapid changes, but 
may as well remain stable over several rotations. These light 
curves have usually been modelled with one simple sinusoid for 
the whole data or for fixed parts of it. In reality, the time span 
of observations required for reliable light curve modelling can- 
not be determined uniquely, because it is certainly not constant. 
Nor does the best model for the light curve remain the same, e.g. 
periodicity can vanish and another periodicity can reappear at 
some later epoch. This means that the selection of the slices of 
data, hereafter called datasets, must be flexible and the complex- 
ity of the model must adapt to the light curve changes. The CPS 
method utilises overlapping datasets, identifies the best model 
for each of these and gives a quantitative estimate for the tempo- 
ral stability of these models. Although the CPS is applied here 
only to simulated and real ground based photometry of active 
stars, this method can be easily adjusted to search for periodicity 
in any type of unevenly spaced data. 

2.1. Datasets and segments 

The input data for the CPS method are observations y, = y(f,) at 
time points f, with observational errors cTi. Before any modelling 
is done, the data must be divided into short datasets (SET), each 
of which is modelled individually. The length of these datasets 
is Ar = f„ - f 1 , which we limit to have the maximum length of 



AT„ 



Po > jAT2, 



(1) 



where Pq is the a priori estimate for the photometric rotation pe- 
riod. The lower limit Ari ensures that most of the datasets con- 
tain enough observations for modelling. Datasets with n < 10 
are not analysed with the CPS. The upper limit AT2 is necessary 
to prevent significant changes in the light curve during a single 
dataset for stars with longer rotation periods. In other words, the 
parameter AT^n-^x is applied to limit the length Ar = f„ - fi of 
the analysed datasets so that they yield a consistent light curve 
model but still have enough observations for reliable modelling. 
In practice, suitable values for ATi and AT2 are 25 d and 200 d 



respectively, when analysing ground based photometry of active 
stars. For different kinds of data, the values of Ari and AT2 must 
be determined separately. 

Unlike in the TSPA, the datasets are allowed to overlap in the 
CPS. This gives a better time resolution. In typical ground based 
observations, we begin a new dataset at the first time point fi of 
each night and include all data points y, = y(ti) within f 1 < f, < 
fi + Ar^ax- A candidate for the next consecutive dataset is chosen 
by moving the window one night forward. We require that the 
observations of two consecutive datasets fulfill the condition 



SETk € SET^+i and SET^+i (t SET, 



(2) 



In other words, both datasets must contain at least one observa- 
tion that does not belong to the other one. If this condition is vi- 
olated, the candidate dataset SET^+i is rejected. In this case, the 
dataset window is moved one night further. If there is much data 
during each night and Po is short, a denser selection of datasets 
might be desirable. In this case, one may even start a new dataset 
at each new observation. This kind of a much denser dataset se- 
lection rule would be ideal for satellite photometry. 

Observations that would be included only in datasets having 
« < 10 are of little use for modelling with the CPS, as for them 
n approaches the number of free parameters of the model. Such 
observations are rejected as temporally isolated data points. 

The data are further divided into longer segments (SEG) 
whenever there is a gap longer than Ar^ax in the observations. 
Such gaps may have contained rejected isolated data points. 
These segments are defined solely for the purpose of identify- 
ing datasets belonging to different observing seasons. 

The benefit of letting the subsets overlap is a better time 
resolution compared to the TSPA, where the datasets were sep- 
arated from each other. As a drawback, the modelling results 
from two overlapping datasets are correlated with each other. 
One can easily eliminate this correlation by comparing only 
the non-overlapping datasets. The alternatives in selecting these 
non-overlapping (hereafter called the independent) datasets are 
numerous. The most simple alternative is to select the indepen- 
dent datasets by beginning from the first dataset of each segment 
and then selecting the next independent dataset with the criterion 
that it does not have any common data points with the preceding 
independent dataset. 

2.2. Modelling of the observations 
The TSPA model 

K 

y{ti) = y{ti,p) = M + 2 [Bt cos {klnftd + Ct sin {klnfu)], (3) 

k=\ 

i.e. a Ki\\ order Fourier series, is also used in the CPS. This 
model has an order K and 2K + 2 free parameters, fi - 
[M, B\,. . . , Bk, Ci , . . . , Ck,/]- In other words, the free param- 
eters are the mean M, the individual cosine and sine amplitudes 
Bj^ and Q and the frequency /. 

In the TSPA, the modelling proceeded through three stages: 
the pilot search, the grid search and the refined search. In 
the CPS, the pilot search, where the initial period estimates 
were identified within a broad frequency interval (Paper I, Sect. 
3.1), is not used. Rather, the grid search (Paper I, Sect. 3.2) 
is performed straight away. In the grid search a dense grid 
of frequencies is tested one at the time by fitting the model 
(Eq. 3) to the observations. When the frequency / is fixed, 
the model becomes linear and the remaining free parameters 



J. Lehtinen et al.: The continuous period search method and its application to the young solar analogue HD 1 16956 



[M, B\, . . . , Bk, C\,. . ., Ck\ of the model have unique solutions. 
The grid search is performed within the period range of 



{\-q)Pa<P<{\+q)Pa, 



(4) 



where Pq is the a priori rotation period estimate of Eq. 1 and 
q regulates the interval of the period search. We use a value of 
q - 0.15. This is usually sufficient, since for many active stars 
there already exist period determinations in the literature and the 
real period changes in photometry can be expected to fall inside 
the ±15% range of Pq- If necessary, this tested period range can 
be expanded. If no Po is available, one can be obtained, for ex- 
ample, with the TSPA. 

The final model parameters p are obtained from the refined 
search (Paper I, Sect. 3.3), which consists of a standard nonlinear 
Marquardt iteration that minimises ;(f^ with weights. 



x\y,~p)^Y,''iel 



(5) 



where w, = ctt^ are the weights and e, = y, - y{ti,p) are the 
residuals. The refined search can find the correct period even if 
it is outside but close to the grid search range (Eq. 4). 

Before modelling, the data is preprocessed by removing out- 
liers. These outliers are determined using a preliminary model 
of order K' . We set this preliminary modelling order equal to 
the highest model order used in the actual modelling of the 
data, K' = K\i^. The preliminary model ^'(0 gives the residuals 
e[ - yi - y'(ti). For each dataset, the observations having residu- 
als larger than three times the standard deviation of all residuals 
e' are removed as outliers. For our data, we found that there are 
typically only a few such outliers among several hundreds of ob- 
servations. 

The problem of determining the order K of the Fourier model 
(Eq. 3) is crucial to the light curve modelling procedure. The 
value of K must be high enough to allow a good fit to the ob- 
served light curve but not too high to result in overfitting to the 
data. In the TSPA, the order K was selected beforehand and the 
same value was used for all datasets. In contrast, the CPS uses 
a Bayesian information criterion to select the best K value sep- 
arately for each dataset. This criterion consists of minimising a 
type of ;^f^ with an additional penalty term for the extra degrees 
of freedom introduced by a higher order K. 

Before testing for the optimal K for each dataset, the up- 
per limit Kiijn for the highest accepted modelling order must be 
selected. Since the smooth light curves of spotted stars usually 
display only one or two minima, not very high values of Kn^n 
are necessary. For more complicated light curve shapes, such as 
those of eclipsing binaries, higher Ku^n values should be used to 
achieve a satisfactory fit. In practice, the value of Kn^n = 2 is 
sufficient for analysing ground based photometry of most spot- 
ted stars. 

The problem of determining K is equivalent to finding the 
number of sinusoids embedded in the data. This number is just 
the order K of the model. A solution for this particular problem 
was presented by Stoica & Selen (2004), who studied the appli- 
cability of different model order selection rules. They argued that 
the best performance is achieved with the Bayesian information 
criterion. It has the property that it always finds the correct K as 
the number of data points n increases to infinity. 

The Bayesian information criterion is given by 



Pbic = 2nlnA(y,p) + (5K + l)lnn. 



(6) 



where the first term gives the logarithmic likelihood with 

A(y,^) - x^iy^P) TIi=i ^i ■ The criterion is used by first mod- 
elling the dataset separately with all models between K = and 
K = Kiiin and then calculating Pbic for each of these models. 
The model with K = (i.e. constant) is obtained simply from 
the weighted mean niy of all data points y, in the dataset, i.e. 
>'(f,) = M = lily. The model for the dataset that minimises Pbic 
has the optimal order K and is chosen as the best model for the 
dataset. 

The function of the two terms in Eq. 6 are the following. 
The first term describes the goodness of the fit to the data. It is 
naturally smaller for higher order models, as they allow a more 
detailed fit. The second term, (5K + l)lnn, is a penalty term 
which grows linearly as a function of K. It balances the smaller 
values of of the first term 2n In A{y,p) at a high K and thus pre- 
vents overfitting. Because of this second term there will be an 
optimal K where Pbic reaches its minimum value and then starts 
to grow again as K increases. 

Most of the free parameters of the model, fi - 
[M, Z?i, . . . , Bk, C\,. . ., Ck, f], are not very useful as such. The 
physically meaningful parameters are the mean magnitude M 
and the total amplitude A of the light curve, the period P - f^^, 
as well as the epochs ?„,!„,! and fmin,2 of the primary and sec- 
ondary minima of the light curve. The mean M and the period 
P = /"' are obtained directly from the free parameters /3 of the 
model, but the values of A, fn,in,i and fmin,2 must be determined 
numerically from the light curve model. Note that only the mod- 
els with K > 2 can provide all of these parameters. For K = I, 
the secondary minimum fmm,2 does not exist and the only param- 
eter that can be obtained for K = Qis the mean M - niy. 

The error estimates and reliabilities of these parameters are 
determined as in the TSPA (Paper I, Sects. 4 and 6.3). This pro- 
cess consists of running a bootstrap with 200 rounds for the 
residuals e and all model parameters. The errors of the param- 
eters are the standard deviations of the 200 bootstrap estimates. 
All these bootstrap distributions are tested against the Gaussian 
hypothesis 

Hg: X represents a random sample drawn from a Gaussian dis- 
tribution, 

where x denotes the residuals or the bootstrap distributions 
of any of the model parameters. If any of these distributions 
fails to satisfy He, all of the model parameters are consid- 
ered unreliable. The Gaussian hypothesis He is tested using the 
Kolmogorov-Smirnov test with preassigned significance level 
y - 0.01 for rejection. Also, if the secondary minimum fmin,2 
is present in less than 95% of the bootstrap samples, it is consid- 
ered unreliable. 

2.3. Time scale of change 

The modelling gives the mean M(t), total amplitude A(t), period 
P(t) and the minimum epochs fmin,i(T) and fmin,2(T) of the light 
curve, where r is the mean of all observing times f, in the current 
dataset. As the shape of the light curve usually evolves with time, 
the model must also evolve. It is of interest to investigate for how 
long an unchanged model for a specific dataset still reasonably 
well fits with the observations of the subsequent datasets. We 
estimate the time scale of change Tq for each dataset i by deter- 
mining for how many of the following datasets l + k the model 
of dataset i still fits the data. 

The dataset i contains the observations y, which yield a 
model yi(ti). Another dataset i + a- at a later epoch contains 
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the data y^+K. Assuming that the model yXti) is valid for both 
datasets, we define the residuals 



ji - yi(ti) 



and 



■y^- 



■ yi(ti+K)- 



(7) 



(8) 



If the light curve has not changed between datasets i and l+k, 
the residuals e^ and e, ^^ should have similar distributions, i.e. they 
should satisfy the hypothesis 

Hk2: £t,K ond e, represent random samples drawn from the same 
distribution. 

Hk2 is tested using the two-sided Kolmogorov-Smirnovtest. The 
preassigned significance level for rejecting Hk2 is y = 0.01. The 
time scale of change is not determined for datasets where the 
model parameters have been found unreliable. This means that 
we already know that the residuals e^ follow a Gaussian distribu- 
tion. Thus the distribution of \ii should also resemble a Gaussian 
distribution, if Hk2 is not rejected. 

The computation of the time scale of change for the model 
of dataset t starts from the dataset l + k, where k - 1 . If the 
sets of residuals e^ and e,^=i for the two datasets pass the test for 
Hk2, the comparison of residuals proceeds to the next dataset, i.e. 
K ^> K+ \, and the same test is applied again. Finally, when this 
test fails, the comparison process is terminated. In this case, the 
time difference between the dataset t and the last dataset passing 
the test, L - K - 1, is taken as the time scale of change for the 
subset L, i.e. 



Tcird 



(9) 



This computation of Tq only proceeds until the end of each seg- 
ment. No datasets from other segments are compared, because 
there is no data confirming or refuting the model during the gap 
between the two segments and this might introduce a bias into 
the Tc estimates. If the computation process hits the end of the 
segment, a special value Tq - -2 is given for the correlation 
time scale. This denotes the case that this particular model de- 
scribes adequately all the datasets following it within the seg- 
ment. 
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Fig. 1. Model order K from analysis of test data simulated with (a) 
model yK=Q{ti) in Eq. 10 and (b) model yK=\iti) in Eq. 10 having am- 
plitude to noise ratio A/N = 10. The number of data points in each 
dataset is presented in panel (c). 

Table 1. Fractions of falsely determined model order K in the analysis 
stable simulated data. 
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K 


A/N = 





A/N = 2 


A/N = 5 


A/N = 10 


A/N = 20 


n<2Q 



1 

2 


0.39 




0.16 
0.17 


0.18 
0.15 


0.12 
0.13 


0.15 
0.16 


n>20 




1 

2 


0.21 




0.10 
0.08 


0.10 
0.09 


0.06 
0.07 


0.07 
0.10 



ti + f 152 - h to the original set of f, and removed the double oc- 
currence of ?i52. This resulted in a set of 303 observation times 
with a complete time span of 377.43 d. 

When simulating the periodic test data, we used a period of 
Po = 7.8 d. This imitated the a priori period estimate for HD 
1 16956. With the chosen observing times and this particular pe- 
riod, our test analyses each yielded a total of 184 datasets with 
^Tjnnx - 25 d, where the number of simulated observations in 
each dataset was in the range 12 < n < 28. 



3. Testing the method 

Before analysing real stellar photometry, we used simulated data 
to determine the performance of the CPS method in certain criti- 
cal situations. This allowed us to get more insight in understand- 
ing some of the results of the analysis. Failing to do this can lead 
to a biased interpretation of the results. 

Since the distribution of the observing times can have a pro- 
found effect to the results of the analysis, we generated our simu- 
lated test data using real observing times of the star HD 1 16956. 
Thus, any spurious effects introduced by the observing times 
alone should also be present in the analysis of both the simu- 
lated data and the real data of HD 116956 (Sect. 4). Because 
other active stars are observed in a similar manner, the results of 
this section are applicable also when analysing them. 

The observing times f, used for generating the simulated test 
data were taken from the first 152 observations of the time series 
of HD 1 16956. These form a complete unbroken observing sea- 
son and are temporally fairly evenly distributed. To get a longer 
time series for the simulations, we appended a duplicated set of 



3.1. Sensitivity to clianges in tine modei order 

As already stated in Sect. 2.2, the Bayesian information crite- 
rion (Eq. 6) should determine the correct modelling order with a 
very high probability when the number of observations is large. 
However, when analysing photometry of active stars, this num- 
ber of observations is often quite limited. Typically we have 
10 < n < 20 in a dataset. This scarcity of data can potentially 
affect the performance of the model order selection. 

To test the performance of the model order selection, we 
analysed three different simulated time series. These consisted 
of one set of constant data and two sets of sinusoidal data each 
with random Gaussian noise 



yK=i(ti) = 0.01 sin(27r/f,) + eN,; 
yK=2(ti) = 0.01 sin(47r/f,) + en,; 



(10) 



Here f ^ - Pq - 7.8 d and eN,, denotes Gaussian noise distribu- 
tion having a zero mean and a variance cr^, i.e. eN,, ~ A^(0, cr^). 
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These three simulation models were chosen so that the correct 
modelling order should he K = 0, K - 1 and K = 2, respec- 
tively. Although both yK=i(td and yK=2(ti) are simple sinusoids, 
the CPS should naturally select K = 2to model the data gener- 
ated with yK=2(ti) when all tested period candidates are close to 
Pq. Then the CPS has to use a higher order model to correctly 
describe the two minima of yK=2{ti) occurring with the period 
Pq. All the analyses were performed with Kn^ = 3. This allowed 
us to examine the possibility that the criterion of Eq. 6 may yield 
a too high K even when analysing simulated data generated with 

yK=2(ti)- 

In the cases of yK=iiti) and yK=2{ti), the test data were sim- 
ulated using four different noise levels where cr-s^ was 0.005, 
0.002, 0.001 or 0.0005. Given the hght curve half amplitude 
A/2 = 0.01, these correspond to amplitude to noise ratios 
A/N = 2, A/N = 5, A/N = 10 and A/N = 20, using the 
definition A/N = Ajlcr^. According to this definition, the data 
simulated with yK=Q{ti) and any ctn has AjN - 0, and thus the 
results for K -Q are independent of the absolute level of noise. 

After analysing the simulated test data, we checked how of- 
ten the CPS arrived at a TT value different from that used in the 
simulation model. Fractions of false K in 15 simulated time se- 
ries (a total of 2760 analysed datasets) are presented in Table 1 . 
They are reported separately for those datasets that had « < 20 
and those that had n > 20. In the case of yK=iiti) and yK=2itd, 
the probability for a false model selection seems to be quite in- 
sensitive to A/N. There is no clear structure visible in the failure 
probabilities at different A/N, nor between K = I and K = 2. 
This seems to indicate that the sensitivity of the selection crite- 
rion does not strongly depend on A/N or the complexity of the 
periodic model. On the other hand, the number of observations in 
the datasets has a substantial influence to its success rate so that 
with more observations the correct model order is detected with 
a much higher probability. Nevertheless, the success rate stays 
between 80% and 95% in both the cases n < 20 and n > 20, 
which can be considered as a satisfactory performance. 

Unfortunately, in the case ofyK=Q(ti), i.e. pure noise, the CPS 
still detected spurious periodicity in 21% of the datasets having 
n > 20 and 39% of those having n < 20. This can cause problems 
if the analysis results are not interpreted critically. Luckily, these 
spurious periodic models can often be identified from their very 
low amplitudes. As a rule, the amplitudes of these models have 
values comparable to the preset noise level os^ and the observed 
standard deviation of the residuals cr^ so that they fulfill A » 
2crN X! 2cr^. This is expected, as the amplitudes are artifacts 
resulting solely from the modelling of noise. 

Examples of the behaviour of K during a test analysis for 
typical simulated data are presented graphically in Fig. 1. The 
analysis of pure noise displays a large scatter of false results, 
although the correct value K = is still clearly the most com- 
monly obtained result (Fig. 1 (a)). The scatter is large especially 
near the mid section of the time series, where the number of 
simulated observations is small, n < 20. The analysis of the data 
simulated with yK=i{ti) and A/N - 10 behaves much better with 
only some overfitting and no false detections of K - models 
(Fig. 1 (b)). 

Clues to why the Bayesian information criterion has trou- 
ble in detecting the correct K for the data simulated using the 
yK=()(ti) model can be seen in Fig. 2. This figure presents sim- 
ulated observations during two arbitrary simulated datasets and 
folded into phase representation using the model period Pq = 7.8 
d. Because of the relatively small number of available data points 
(n - 19), the case of pure noise, i.e. underlying K = 0, still 
vaguely resembles some periodic form (Fig. 2 (a)). The proba- 




Fig. 2. Examples of single 25 d long datasets simulated with (a) model 
>'A:=o(fi) in Eq. 10 and (b) model >'r:=i(/,) in Eq. 10 having amplitude to 
noise ratio A/N =10. The phases of these simulated observations are 
calculated using the period of the simulation model, Po = 7.8 d. 




200 
time [d] 

Fig. 3. Primary (squares) and secondary (triangles) minimum phases 
from analysis of test data simulated with >'spoi(f,7/'o) of Eq. 13. Filled 
symbols denote reliable and open symbols unreliable phase estimates. 
The two straight lines denote the correct phases where the simulated 
spots should have been detected. 



bility of chance resemblance of periodic data decreases only for 
a larger number of data. In the case of an underlying K - I sim- 
ulation model (Fig. 2 (b)) it is quite clear that a single sinusoid 
is the correct model for the data. 

We conclude that the CPS detects real periodicity even from 
low amplitude data of a high noise level. If there is no real pe- 
riodicity, spurious periodicity is still detected in about a quarter 
of the analysed datasets. These spurious period detections can, 
however, be identified because their observed amplitude to noise 
ratio is near unity. 



(A/N)obs = :^ * 1- 

In conclusion, the results need always to be interpreted with spe- 
cial care when the dataset contains few observations and the 
model amphtude A is low. 

3.2. Sensitivity to two dose minima 

When analysing the distribution of starspots on the surface of a 
star, one often simply identifies the light curve minima as direct 
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Fig. 4. Two single spot model light curves according to >'spot,i (<^) plotted 
around phases (j)o = 0.15 and (po = -0.15 (dashed lines) and the com- 
bined light curve (solid line). The phase separation A<p is indicated with 
the two headed arrow. 



markers of the spot longitudes. In the case of a single spot or two 
clearly separated ones, this assumption, of course, holds and the 
light curve minimum phases give good approximations of the 
spot longitudes. 

However, if two spots are longitudinally close enough to 
each other, the observable light curve minima can be signifi- 
cantly shifted. As the two minima produced by the two spots 
move closer to each other, their tails start to merge. At first, the 
minima stay separated but are shifted somewhat towards each 
other. At phase separations below some critical value A^cdt, the 
minima finally merge and produce one single broad minimum. 
The observed phase of this single minimum lies between the 
phases of the two underlying spots. 

The tendency of the light curve minima, produced by two 
individual spots, to merge depends on their width. Note, that no 
minimum can be very narrow, even in stars with an inclination 
close to 90°. For example, a spot located at the stellar equator 
is typically visible over half of the stellar rotation. During this 
time it modulates the visible brightness of the star. Thus the total 
width of the minimum spans about one half of the total rota- 
tion period. If the spot lies closer to the visible pole or the stel- 
lar inclination is small, it stays at the visible stellar disk for a 
longer period and causes an even broader minimum. In the ex- 
treme case, a spot can stay at the visible hemisphere at all times. 
In this case it causes a very broad minimum induced only by 
its variable projected area and limb darkening. Only a spot be- 
low the equator, closer to the unseen pole, may cause a narrow 
minimum, as it stays visible for less than half of the rotation pe- 
riod. However, all minima significantly narrower than half of the 
rotation can only be produced by spots that appear on the visi- 
ble hemisphere briefly and never come very far from the limb of 
the star. Because of limb darkening and a small projected spot 
area, these minima are very shallow and are therefore strongly 
affected by noise. 

To get an estimate of the critical phase separation A0crit, con- 
sider a model light curve of the form 



fO, 

y^potA<f>) = \ 5[1 -H cos (4710)] 

0, 



(/) < -0.25 
-0.25 <<p< 0.25 
d> > 0.25 



(11) 



in phase space within the interval -0.5 < (p < 0.5. This model 
approximates reasonably well the shape of a light curve pro- 
duced by a single spot at the stellar equator at phase (p - Q. 
Putting two equally strong spots at phases (po and -(po results in 



a merged light curve. The merged light curve within the interval 

(pQ - 0.25 < (p < 0.25 - (pQ has the form 



Note that this is only a linear approximation, as actual photom- 
etry is measured in magnitudes which can not be added in this 
manner themselves. In the case of starspot induced light vari- 
ability where the amplitude of the light curve is in the range of 
0.01 - 0. 1 mag, this linear approximation can, however, be done. 
The combination of two 3'spot,i(0) curves at the phases 0o = 0.15 
and 00 - -0.15 is visualised in Fig. 4. The two underlying 
curves are drawn with dashed lines and the merged curve with a 
solid line. The phase separation of the two components is indi- 
cated with the two headed arrow. 

The two spots can be identified from the light curve if they 
cause two separate light curve minima. In other words, there 
should be a secondary maximum at phase (p - Q between the 
two minima. This corresponds to the condition 



^sUaCO) = -87r2[cos(47r(0o)) + cos (47r(-0o))] < 

^ cos (47r0o) > 0. 

For < (pQ < 0.25, this corresponds to the phase separation of 
the two spots being 



^4> = 200 > 0.25, 



(12) 



i.e. no spots closer to each other than A0crit = 0.25 should cause 
two separate light minima. Less than optimal performance of the 
analysis algorithm may yield an even larger value of A0crit- 

To determine the phase resolution of the CPS, we tested a 
two spot simulation model based on two adaptions of Eq. 11, 



yspoti- -yspouh^ +3'spouU,,. 



^p 



+ 0.1 + EN,,. (13) 



The model describes two spots rotating with periods Pq and 
Po + AP and having an initial phase separation of A0 - 0.1. 
The difference in period, AP = -0.001327Po, was chosen so 
that at the end of the 377.43 d long time series the phase differ- 
ence between the two spots would have increased from A0 = 0. 1 
to 0.6. The random noise used in simulating this time series was 
fN,; ~ A^(0,0.00l2), i.e. A/N = 10. 

The resulting minimum phases from the CPS analysis of this 
simulated data are shown in Fig. 3. Superimposed on the mini- 
mum phases detected with the CPS are two straight lines indicat- 
ing the correct simulated phases where the two spots should have 
been found. At the beginning of the time series, the two under- 
lying minima are inseparable and produce one common merged 
minimum with a phase between the two correct phases, as ex- 
pected. Towards the end of the time series, the two minima both 
become detectable and have estimated phases very close to their 
correct values. However, near the time when the two minima first 
become separable their estimated phase differences are system- 
atically smaller than the correct simulated values. 

The two minima are first detected separately at time t - 
167.99 d at which time their phase separation given by the CPS 
is A0 = 0.27. The simulated correct phase separation at that mo- 
ment is A0 = 0.32. Both of these two values exceed the limit 
A0crit = 0.25 derived in Eq. 12. This indicates that the ability of 
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Fig. 5. Period estimates from the analysis of test data simulated with 
yK=i{ti) having A/N = 10 (see Eq. 10). 

Table 2. The spurious changes of the period P in units of Z (Eq. 14) for 
test data simulated withyn^iitj) (Eq. 10) and different values of A/N. 
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the CPS to separate individual minima is not perfect. We con- 
clude that spots with a phase separation smaller than about one 
third of a rotation are not expected to be separable with the CPS. 
The failure to detect two close by spots separately from a 
light curve means that any time when two spots are detected 
from a light curve, they already have a large longitudinal sep- 
aration. In other words, our impression of the distribution of the 
spots is biased towards spots lying nearly at opposite sides of 
the star. This can lead to detections of spurious active longitudes 
even in the case when the underlying physical spot distribution 
has a completely different geometry. 

3.3. Stability of tine metiiod 

The results of the analysis should be reliable regardless of ran- 
dom effects caused by noise and uneven distribution of the ob- 
serving times. Analysis of test data simulated with an unchang- 
ing underlying model should thus give stable parameter esti- 
mates throughout the whole time series. 

To test the stability of the CPS, we analysed sinu- 
soidal test data simulated with yK=i{ti) (Eq. 10). The stan- 
dard deviation of the noise was selected from crj^ e 
[0.0005, 0.0007, 0.001, 0.002, 0.003, 0.005] corresponding to 
amplitude to noise ratios A/N = 20, A/N = 14.3, A/N = 10, 
A/N = 5, A/N = 3.3 and A/N = 2, respectively. The analy- 
sis yielded generally rather stable values for the simulated light 
curve mean M and amplitude A and the single minimum phase 
0min,i- However, the period estimate P showed quite significant 
variations, as can be seen in Fig. 5. This type of spurious vari- 
ations must be taken into account, since the period changes of 
spotted stars are commonly used to measure the strength of their 
surface differential rotation (Jetsu 1993). 

To quantify the effect of the spurious period variations, we 
calculate the parameter 
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Fig. 6. The differential photometry of (a) HD 1 16956 and (b) the com- 
parison star HD 1 14446. Each observing season is labeled with its seg- 
ment number SEG=1,2,...,12. The filled triangles connected with con- 
tinuous lines denote seasonal mean differential magnitudes. 



from the period estimates of the simulated data. This parame- 
ter measures the variability of P within its weighted ±3cr lim- 
its. The definitions of the parameters are w, = cr^^, P„ = 

2:w,P,[2:>v,r' and APw = VSw^K^^^^A^/ Vr^ (Jetsu 
1993). The Z values calculated from the independent datasets 
for all simulated models are given in Table 2. When using the 
same parameter Z to characterise the differential rotation of a 
star, we may call these values spurious differential rotation Zspu. 

What is striking in the results of Table 2, are the large val- 
ues the spurious differential rotation Zjpu, especially with lower 
A/N. This severely affects the detectability of weak differential 
rotation. For example, real differential rotation causing variabil- 
ity of only 2% in P would be hard to distinguish from the spu- 
rious differential rotation even from data of good quality with 
A/N - 20. The values of Zspu are roughly inversily proportional 
to A/N. Even very large 15% variations of P would be indistin- 
guishable from the spurious differential rotation with A/N = 2. 

By considering the period variations caused by physical and 
spurious differential rotation, i.e. Zphys and Zjpu, to be indepen- 
dent effects, we can estimate the contribution of the real physical 
differential rotation. In this case, the observed value of Z^ is the 
sum of the squares of the two components, Z^ 
yielding 



■^phys ^ -^spu' 



ryl _ 72 _ 72 

■^phys ~ -^ -^spu- 



(15) 



(14) 



The value of Zspu in this equation can be interpolated for any 
A/N based on the values in Table 2. In theory, Eq. 15 gives an ex- 
act instantaneous value for Zphys under the assumption that Zphys 
and Zspu are uncorrelated. Note, however, that both the ampli- 
tude A of the light curve and the observational accuracy a-,, are 
likely to change during a time series. Thus it is impossible to get 
a unique estimate for Zspu- Therefore, the relation of Eq. 15 only 
gives a rough guideline for estimating the effect of the spurious 
period variations. 



4. Analysis Of HD 116956 

In this section we present results of the CPS analysis of long- 
term photometry of the young solar analogue HD 116956. The 
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analysis was done using the a priori period estimate Pq = 7.80 d 
based on the preliminary results of Gaidos et al. (2000). The 
upper limit for the modelling order was set to be Knja - 2. 

To further justify the choice of Pq, we can estimate the lower 
limit for the stellar radius, Rsini - Pvsmi/2n = 0.86/?©, by 
using P = Pq and vsin/ = 5.6 km s"' as given in Sect. 1. This 
limit is slightly smaller than the radius of the Sun, as expected for 
a star with a somewhat later spectral type. Another estimate for 
the stellar radius can be obtained using the Barnes-Evans relation 
(Lacy 1977) 



log (R/Re) = 7.4724 - 0.2yo - 2Fv + log d, 

where Fy = 3.977 - 0.429(y - R)o and [d] = pc. Using 
Vo = 7.286, {V - R)q = 0.5 and <i = 21 .9 pc from the Hipparcos 
(Ferryman et al. 1997) and USNO-B (Monet et al. 2003) cata- 
logues and neglecting the interstellar reddening because of the 
small distance to the star, we get R = (0.66 ± Q.06)Rq. This value 
is slightly smaller than the lower limit calculated using the ro- 
tation period and velocity. Nevertheless, the two results for the 
radius are of the same order and Pq is a reasonable estimate for 
the rotation period. 

Our photometry' of HD 116956 was obtained over 12 
consecutive observing seasons between HJD = 2451172 (24 
December 1998) and HJD = 2455342 (25 May 2010) with the 
T3 0.4 m automatic photoelectric telescope (APT) at Fairborn 
Observatory in Arizona. The APT performs differential photom- 
etry in the Johnson V andBpassbands. Abrief description of the 
operation of the APT and the reduction of the data can be found 
in Fekel & Henry (2005) and references therein. A total of 1408 
differential V band observations were acquired with HD 1 14446 
as the comparison star This data is shown in Fig. 6, which also 
shows the season numbers that coincide with the segment divi- 
sion of the CPS analysis. The seasonal means are denoted with 
filled triangles which have been connected with continuous lines. 

The external precision of the observations taken with the T3 
telescope is 0.004 - 0.005 mag, as shown by simultaneous mea- 
surements of the comparison star HD 1 14446 against a check 
star HD 119992 (see Fig. 6 (b)). The difference between the 
comparison and check star showed no seasonal changes during 
the whole observing period, i.e. there are no prominent system- 
atic errors in the data. This is also supported by the ;^f^-test for 
the n = 1273 comparison star minus check star differential mag- 
nitudes. Using e; = y; - my and cTi = 0.005 gives ;t'^ = 1166.6 
(Eq. 5). This corresponds to the critical level Q = 0.9829, where 
Q states the probability of ;^f^ reaching the observed value un- 
der the null hypothesis that these differences remain constant. 
There is no need to reject this constant brightness hypothesis. In 
contrast, the n = 1408 HD 116956 minus comparison star dif- 
ferential magnitudes give x^ - 12002.9 which corresponds to 
Q <sc 10"^", i.e. these observations certainly exhibit variability. 
For a more detailed description of the data acquisition see Henry 
(1999). 



4.1. Graphical presentation of ttie results 

The CPS divided the data into twelve segments, each constitut- 
ing of a complete observing season. The observing seasons typ- 
ically lasted from the beginning of December to the beginning 
of July. The first season of our data has also been analysed by 



' The V band photometry used in this paper is published electroni- 
cally at the CDS. 



Gaidos et al. (2000). The analysis of this first segment (SEG = 1) 
is presented graphically in Fig. 7, where the units are V magni- 
tudes for (Tf , M and A and days for Tc and P. All the rest of the 
parameters are dimensionless. This figure contains the subplots, 
where the parameters are plotted as functions of t: 

(a) Standard deviation of residuals cr^(T) 

(b) Modelling order K(t) (squares, units on the left y-axis) and 
number of observations per dataset n (crosses, units on the 
right y-axis) 

(c) Mean differential V-magnitude M(t) (HD 1 16956 minus HD 
1 14446 

(d) Time scale of change rc(T) 

(e) AmpHtudeA(T) 

(f) Period P(t) 

(g) Primary (squares) and secondary (triangles) minimum 

phases 0min,l(T) and 0min,2(T) 

(h) M(t) versus P(t) 
(i) A(t) versus P(t) 
(j) M(t) versus A(t). 

In the subplots (a), (c) and (e)-(g), the reliable parameter es- 
timates are indicated by filled symbols and unreliable ones by 
open symbols. The reliability of the parameter estimates is de- 
termined as described in Sect. 2.2. In the case of P, also the 
level of the a priori period estimate Po is shown (horizontal line). 
This fits reasonably well to the detected P values since the same 
data of SEG = 1 was used by Gaidos et al. (2000) to determine 
P(). The minimum phases are calculated using the median pe- 
riod Pmed of the segment. In subplot (d), the upward pointing 
arrows indicate that the computation of the time scale of change 
has reached the end of the segment and they correspond to the 
values Tc = -2 (see App. A). In the correlation plots (h)-(j), 
the error bars have been drawn only for the parameter estimates 
of independent datasets. The linear Pearson correlation coeffi- 
cients ro for the independent datasets, as well as the probabili- 
ties P(\r\ > ro), are given. Finally, the values of the a priori period 
estimate Pq, the median period Pmed, the limiting modelling or- 
der Kiiin and the maximum length the a dataset Ar^ax are given 
above the plot. 

The light curve fits for the independent datasets are pre- 
sented in Fig. 8. Because the CPS models each of the datasets 
with a different period P(t), it is in principle impossible to rep- 
resent these light curves with an ephemeris based on one con- 
stant period. Each dataset was therefore first modelled with the 
phase (pi = FRAC[(f - fmin,i(T))/P(T)], i.e. a different period 
value was used for each dataset. Then the phases (p-^n of the pri- 
mary minimum epochs tj^inj were computed with the constant 
period ephemeris HJDai = 245 1176.94+7.8416£' (see Sect. 4.3). 
Finally, the data and the light curve models were displayed as a 
function of the phase (p - (pi + (pai,\- We chose this particular 
definition for the phase, because it allows an easier comparison 
between Figs. 8 and 10. 

4.2. Long term variations 

Our Fig. 9 displays the long term variation of the light curve 
parameters M, A and P. With the exception of the rotation period 
P, these seem to display some regular behaviour. This may be a 
sign of activity cycles in the star. The variations of P seem to be 
more or less random. 

The long term variations are most striking in the mean 
brightness M. During 2004-2005 it underwent a dip of ~ 0.02 
mag. A similar dip occurred later in 2008. Also at the beginning 
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Fig. 7. The CPS analysis of SEG = 1 of the photometry of HD 116956. Descriptions of the subplots are given in Sect. 4.1. 



of the observations in 1999, the values of M were dimmer than 
the average suggesting another similar dip. The same variability 
can be seen in the raw V data (Fig. 6). In addition to the minima 
before 1999 and in 2005 and 2008, there seems to have been one 
during 2001-2002. This one cannot be seen in the M variations. 
It is possible that signs of this particular dip are lost in the gap 
between the observing seasons of 2001 and 2002. Put together, 
these timings might suggest an activity cycle of roughly 3 years 
in length. 

The variations of the mean brightness of the star have the 
most straight forward interpretation as an activity proxy. As the 
spottedness of the star increases, its mean effective temperature 
decreases and a lower brightness is observed. Another proxy is 
provided by the light curve amplitude A. It measures more or less 
the nonaxisymmetry of the spot distribution. There may or may 
not be a correlation between this and the total spottedness and 
likewise between the variations of A and M. As it happens, the 
connection between the long term variations of M and A is not 
immediately clear. There seem to be variations in A with about 
the same time scale as in M. However, these variations are more 
erratic than those of M. 



To investigate somewhat further, whether the star has any ac- 
tivity cycles, we performed the CPS analysis for the independent 
M and A estimates, using cr'^ and cr^^ respectively as weights. 
Our a priori cycle period estimate was Pq = 3.0 yrs. The result- 



ing cycle periods were Pm = 3.26±0.04yrs andPA = 3.09±0.13 
yrs. However, these models had;^f^ = 1206.2 » «ind = 72 and 
x\ - 847.1 » Hind - 72, where nind is the total number of 
the independent datasets. Any reasonable model should satisfy 
X^ ~ "ind, i-C- these two cycles are simply not statistically signif- 
icant. Thus, the data does not seem to support the interpretation 
of the seasonal M and A variations as signs of activity cycles. 

There may be short time scale correlations between M, A 
and P. For SEG = 1, the linear Pearson correlation coefficients 
are ro(P,M) = -0.649, ro(P,A) = 0.324 and ro(A,M) = -0.179. 
Because of the small number of independent datasets («ind = 
7) none of these are significant. The same applies for all other 
segments as well. Other than simply linear correlations are also 
possible, though even more difficult to detect. The presence of 
such correlations is suggested by paths of parameter estimate 
pairs in Figs. 7 (h)-(j). 



4.3. Active longitudes 

The long term variation of ffie primary and secondary light curve 
minimum phases is displayed in Fig. 10 with the same notation 
as in Fig. 9. The minimum phases (pmm,i and (pmm,2 are calculated 
from the minimum epochs fmin.i and fmin,2 using the period P^i = 
7.8416 ± 0.0011 d. This period was detected when the Kuiper 
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Fig. 8. Light curves of independent datasets: Each dataset was first modelled with the phases 0i = FRAC[(/- rmin,i(T))//'(T)]. The phases 0i,i.i of 
the primary minimum epochs fmin.iCi") were then computed with the constant period the ephemeris HJDai = 2451176.94 + 7.8416ii (see Fig. 10). 
Finally, the data and the light curves were plotted as a function of the phase (p = (pi+ 0j,i i . 
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Fig. 9. Evolution of M, A and P during the complete time series. 
Parameter estimates from independent datasets are denoted as squares 
with error bars and from all other datasets as small crosses. M and A are 
given in V magnitudes and P in days. 



test was applied to the n - 90 reliable primary and secondary 
minimum epoch estimates of all independent datasets. 

The formulation of this nonweighted K-method was the 
same as in Jetsu & Pelt (1996). It tests the null hypothesis (Ho) 
that the phases of the epochs of the primary and secondary min- 
ima are evenly distributed between and 1 . Under Ho, one can 



determine the probability that the K-method test statistic reaches 
any particular value within the chosen tested period interval. 
This probability is called the critical level Qk and if it is very 
small, the null hypothesis must be rejected. The best detected 
period has the smallest Qk- 

The period interval that was tested with the K-method was 
between Pmi„ = 0.85Po = 6.63 d and P^ax = l.lSPo = 8.97 d. 
The critical level of the best 7.8416 d periodicity was extremely 
significant, Qk = 8.7 ■ 10"". It exceeds all critical level esti- 
mates given in Table 2 by Jetsu (1996), where active longitudes 
were detected in A And, a Gem, II Peg and V711 Tau. Such 
an extreme value of Q^ confirms the presence of two long-lived 
active longitudes in HD 1 16956. 

A striking feature of the minimum phases is that they are 
confined to two active longitudes with a phase separation of 
A(f> =s 0.5. Comparing Figs. 10 and 7 (g) shows that this phe- 
nomenon is present in both long and short term evolution of 
the light curve. Considering the results of the analysis in Sect. 
3.2, the existence of two confined areas of light curve minimum 
phases could be attributed to a selection effect. However, both 
the long term stability and the fact that there is very little short 
term random changes in the minimum phases show that the ob- 
served active longitudes are indeed likely to be real. If they were 
created by a selection effect, they would not remain stable on 
time scales longer than the length of single datasets Ar^ax- 

The existence of a stable period with which the active longi- 
tudes rotate indicates that there is some coherent magnetic struc- 
ture inside the star rotating with the period Pal- Such a structure 
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ephemeris HJD,, = 2451176.94 + P.i^, where Fj = 7.4816 days. 
Estimates from independent datasets are denoted as squares (primary 
minima) and triangles (secondary minima) with error bars. Estimates 
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could be generated by a nonaxisymmetric dynamo mode (Krause 
&Radlerl980,p.271). 

In addition to the long term stability, the active longitudes 
seem to have experienced periods of migration. This is evident 
especially from the seasonal primary minimum phases. Between 
the observing seasons of 2004 and 2005, there was a phase jump 
of A(p ~ 0.3. During the next years, the primary minimum mi- 
grated to the same phase where it had resided before the phase 
jump. Quite remarkably, this transient occurred at the same time 
as the minimum of M at 2005. It is not clear whether these two 
phenomena are related to each other. Another short phase jump 
seems to have occurred at 2009 and possibly one also at 2000. 
In both cases, the primary and secondary minima seemed to ex- 
perience the same phase shift with their mutual phase separation 
staying more or less constant. 

Another remarkable property in the distribution of the mini- 
mum phases is that for most of the time the primary minima are 
confined to one active longitude and the secondary minima to the 
opposite one. This is evident in Fig. 10 where the primary min- 
ima present in the independent datasets are denoted as squares 
and the secondary minima as triangles. Only at the beginning of 
the time series are there exceptions to this. Two of the indepen- 
dent datasets had their primary minima near the phases where 
the following independent datasets had their secondary minima 
and one had its secondary minimum near the phases where the 
other independent datasets had their primary minima. The course 
of events can more closely be examined in Fig. 7 (g). At the 
beginning of the first segment, the primary light curve minima 
were located near - 0.7. As the segment continued, there was 
a migration of observed primary minimum phases to (f> = 0.2, 
where a new active longitude formed. Roughly 70 d after the 
start of the segment, there was a shift in the order of the min- 
ima so that the secondary minimum grew deeper and became 
the new primary minimum. Towards the end of the segment, the 
old primary minimum gradually faded away. This change from a 
primary minimum to a secondary and vice versa if a direct proof 
of the occurrence of a flip-flop event in HD 1 16956 (Jetsu et al. 
1993). 

Another closer look at the flip-flop in segment 1, that is ob- 
servations between HJD = 2451172 and HJD = 2451361, is 



Fig. 11. The light curve of HD 116956 during SEG = 1 folded with 
median period P^ed = 7.8896 d and divided into three parts. Panel (a) 
includes observations from the first 29 nights of the segment, panel (b) 
observations from the subsequent 56 nights and panel (c) observations 
from the last 104 nights. The scale of the Icr accuracy of the data (0.005 
mag) is indicated with a separate error bar at the right upper corner of 
the panels. 



provided by Fig. 11 . It displays the observations folded with the 
median period Pn,ed = 7.8896 d of the segment. Because of pe- 
riod variations during the segment, the resulting light curves are 
blurred. Nevertheless, the general variations of the light curve 
shape remain clearly visible. The observations in Fig. 1 1 are di- 
vided into three parts, panel (a) includes observations from the 
first 29 nights of the segment, panel (b) observations from the 
subsequent 56 nights and panel (c) observations from the last 
104 nights. During the first and second parts of the segment, the 
light curve displayed two minima near - 0.2 and (p = 0.8 and 
its total amplitude remained relatively low. A distinguishing fea- 
ture between the two parts is the relation between the depths of 
the two minima. During the first part, the minimum near (p = 0.2 
is distinctly the deeper of the two. During the second part, no 
clear differences between the depths are visible. Quite remark- 
ably, there was an abrupt change in the light curve shape be- 
tween the second and third parts. The minimum near = 0.8 
was slightly shifted and the minimum near = 0.2 disappeared. 
Also the light curve amplitude increased at the same time. As the 
observations continued unbroken during the changes, we must 
conclude that the changes were relatively fast. At the same time 
with the changes in the light curve shape there was a jump in 
M (Fig. 7 (c)), which may have been connected to this flip-flop 
event. 

The phase migrations at 2005 and 2009 may not be inter- 
preted as flip-flops in the sense of Jetsu et al. (1993), as there 
were no switches of activity between the two active longitudes. 
Rather, both the primary and secondary minima remained at 
their own active longitudes, as well as preserved their phase sep- 
aration. 
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4.4. Differential rotation 

If the stellar photosphere rotates differentially, we can expect to 
observe different rotation periods for spots at different latitudes. 
At any particular time, only one value of P can be observed. But 
as the spot distribution evolves, so should the observed rotation 
period. The observed period values correspond to the total con- 
tribution of all the spots on the star weighted by their contrast to 
the unspotted surface and their visibility. 

As already described in Sect. 3.3, we can measure the 
strength of the surface differential rotation using the ±3cr lim- 
its of P variations (Eq. 14). Using only period estimates from 
the independent datasets, we get the weighted mean period 
Pw ± AP„ = 7.8288 ± 0.1423 d and the strength of the varia- 
tions Z= 0.11 = 11%. 

Combining the observational accuracy of 0.005 mag and a 
typical light curve amplitude A - 0.05 mag gives a signal to 
amplitude ratio A/N = 5. Using the results of the analysis in 
Sect. 3.3, this A/N value would introduce spurious period vari- 
ations of Zspu = 0.06 = 6%. We can investigate the contribu- 
tion of this effect by calculating a rough estimate for the physi- 
cal component of the period variations using Eq. 15. This gives 
•Zphys ~ 0.09 = 9%, which is still larger than Zspu, although being 
of the same order of magnitude. Thus there remains a possibility 
that even these period variations are not physical and no differen- 
tial rotation is observed. The possibility of a constant period in 
the observational data can, however, be rejected by computing 
the x^ value (Eq. 5) of the tiind - 72 independent period esti- 
mates against the hypothesis P - P^ = constant. The residuals 
ei - Pi - Pw and the errors crp. give x^ - 282, which corre- 
sponds to a critical level Q <s 10"^", i.e. the hypothesis of a 
constant period must be rejected. 

Assuming that the period variations are indeed caused by dif- 
ferential rotation, we can try to interpret the Z value in somewhat 
more depth. If we make the typical assumption of solar like dif- 
ferential rotation profile P{b) - Peq/(1 - ksin^b), where b is 
latitude and Pgq the equatorial rotation period, we can relate the 
calculated Z value to the differential rotation coefficient k with 
the scaling law \k\ a; Z/h, where h - sin^ batax - sin^ femin and 
foniin and /jmax are the minimum and maximum latitudes of spot 
formation (Jetsu et al. 2000). In order to do this we must, how- 
ever, know the latitudinal extent \bmm, ^max] of the spot activity. 
If the spots are restricted to a narrow latitude interval, we get 
observational data only from the rotation periods within that re- 
gion. The actual value of A: is thus larger than what Z suggests. If 
the spot activity spans all the way from the equator to the poles, 
the scaling coefficient h approaches unity and Z » \k\. However, 
spots near the poles contribute less, or not at all, to the rotational 
modulation of brightness. This means that the suitable values of 
h are always somewhat less than unity. 

In the case of HD 1 16956 we have no knowledge of the lati- 
tude extent of the spot activity. We can, however, apply different 
scaling factors. In the case of total range of spot activity from 
fomin = 0° to /7max = 90°, wc havc /i = 1 and |/t| = 0. 1 1 . This is 
roughly half of the solar value of A: = 0.20. But for a solar like 
distribution of spot activity from b^y^ - 0° to /j^ax = 30°, we 
get h ss 0.25 and \k\ - 0.44, over twice the solar value. Such a 
large value renders the assumption of solar like spot distribution 
highly dubious. It would be more likely that the spot activity on 
HD 1 16956 is distributed over a larger latitudinal range than on 
the Sun and the value of fc be closer to A: = 0.1 1. This would cor- 
respond to diflFerential rotation rate AQ = kQ. - InkjP - 0.088 
rad d ' between the equator and the poles. 



Henry et al. (1995) investigated the relation between the ob- 
served values of k and the stellar rotation periods P^^t and arrived 
at the relation 



logfc= -2.12 + 0.76 log P„t-0.57F, 



(16) 



where Piot is in days. The Roche lobe filling factor F - 
^star/^Roche rcduccs lo F - Q for single stars. Using the weighted 
mean period P^ - 7.8288 d we get a prediction k - 0.036 for 
HD 1 16956. This is a low value, about one third of the observed 
lower limit k > 0.11. Thus it seems that, both the lower and 
upper limits k - 0.11 and k - 0.44 reside in or near the scat- 
tered region of differential rotation estimates in Fig. 28 of Henry 
et al. (1995), and therefore HD 1 16956 would appear to undergo 
much stronger differential rotation than expected. 

Another comparison was done using the results presented by 
Collier Cameron (2007). Their relation between diflFerential ro- 
tation rate (AQ) and the eflFective temperature of the star (Tetf) 



AQ = 0.053 



\5130/ 



(17) 



where Teff is in kelvins and AQ in radians per day. This predicts 
AQ = 0.057 rad d ' for HD 116956. This is again smaller than 
the observed value, although closer to it than what was predicted 
earlier by Eq. 16. 



4.5. Time scale of change 

The parameter Tc is an estimate of the typical time scale in 
which the spot distribution changes. It is an important param- 
eter when investigating the short term dynamics of the spot ac- 
tivity. Unfortunately, interpreting the values of Tc is not simple. 
As can be seen in Fig. 7 (d), the Tq values evolve quite rapidly. 
Moreover, there are a large number of datasets where the model 
has been applicable right to the end of the segment (denoted with 
arrows), i.e. no Tc estimate has been obtained. 

The seemingly random evolution of Tq stems from many 
causes. First of all, the complexity of the model aflFects the time 
that this same model stays applicable to future observations. 
Simple models with low K tend to have longer Tc than more 
complex ones. This is apparent from the very first few datasets 
of SEG = 1, as they have K = 1 (Fig. 7 (b)) and Tc > 160 d 
(Fig. 7 (d)) reaching the end of the segment. 

An incomplete phase coverage of the light curve can also 
cause a too low value of Tc- If there is a phase gap in the light 
curve with no observations, there may be a situation where the 
step to the next dataset introduces new observations to the phase 
gap thus completing the phase coverage. In this case, the contri- 
bution of new observations can significantly alter the light curve 
so that the model of the previous dataset is no longer applicable 
and Tc gets a very low value. 

Even if there are no significant phase gaps in the observa- 
tions, an increasing amount of observational data can aflFect the 
value of Tc. This is because with more data the light curve can 
be determined more accurately and there is a higher probability 
of choosing a more complex model. This leads to smaller val- 
ues for Tc. Through a concrete example, it can be seen how a 
larger number of observations and a better observational accu- 
racy decrease the values of Tc. These effects can be nicely seen 
in the observations of e Eri by Croll et al. (2006). They analysed 
continuous photometry of the star taken with the MOST satel- 
lite during three rotation periods. Their light curve is both well 
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sampled and accurate and it turns out to have a different shape 
during each of the three rotations, as can be seen in their Fig. 2. 

Because of the complex behaviour of Tc, it is best to look 
at the long-term mean time scale of change rather than studying 
the short term changes. This long-term mean for HD 116956 
is Tc = 44.1 d. This value is nearly two times as large as the 
length of the datasets Ar^ax = 25 d. Thus the light curve usually 
remains unchanged during one typical dataset and the chosen 
value of ATnjax is reasonable. 

The timescale of light curve change can be compared to 
the convective turnover time t^, as mixing of the convection 
zone may alter the spot distribution which is then observed as 
a change in the light curve. Kim & Demarque (1996) performed 
theoretical modelling to find the values of Tc in late type stars. 
Ossendrijver (1997) published an interpolation formula for their 
results by fitting a cubic polynomial to the theoretically calcu- 
lated values, 



Tc = -68.3 -t- 224.8x - 177.2x2 ^ 57.0x^ 



(18) 



where x - B - V. This formula gives t^ in days. By using the 
Hipparcos value of B - V = 0.823 for HD 116956 (Ferryman 
et al. 1997), we get Tc = 28.5 d. Even if the two values are of 
the same order of magnitude, it is not clear whether there is a 
meaningful link between t^ and Tq- 

5. Conclusions 

We have formulated a new Continuous Period Search (CPS) 
method by improving the earlier TSPA method described in 
Paper I. Our goal was to develop better tools for the analysis 
of photometry of active stars. The three new features are: 

(1) A sliding window for selecting overlapping datasets. 

(2) A criterion for selecting the correct order for the model. 

(3) A time scale for the change Tc of the model. 

Using the sliding window significantly improves the time reso- 
lution of the analysis. The CPS determines the light curve pa- 
rameters with a time resolution much shorter than the length of 
individual datasets. In contrast, the TSPA could only achieve a 
time resolution equal to the dataset length. A similar sliding win- 
dow was applied for a first order harmonic model in Berdyugina 
et al. (2002, their Eq. 3). Our model order selection criterion 
allows the selection of a model with suitable complexity for de- 
scribing the data in each individual dataset. This even allows the 
alternative of a constant light curve with no periodicity. Since the 
CPS can use Fourier series of any arbitrary order as a model, it 
can in principle be applied to any type of periodic data. Finally, 
the computed Tc parameter measures the time scale of change 
of the light curve. 

In order to characterise the performance of the CPS method, 
we tested it with simulated data. We found imperfect perfor- 
mance in three critical situations. 

Firstly, the model order selection criterion does not give the 
correct order K in all of the analysed datasets. This is caused by 
the limited number of observations. Typically, when analysing 
photometry of active stars, each dataset contains 10 < n < 30 
observations. The success rate in finding the correct modelling 
order depends on n. In the case of analysing periodic data, the 
success rate is between 80% and 90% for datasets with n < 20 
and larger than 90% for datasets with n > 20. In the case of con- 
stant noisy data with no periodicity, this success rate is much 
smaller, and spurious periodicity is still detected in 40% and 



20% of the datasets having n < 20 and n > 20, respectively. 
This reduces our ability to distinguish between time intervals 
of periodic variation and constant brightness in the light curve. 
However, this type of spurious periodicity can often be identified 
from the low observed amplitude to noise ratio of the models. 

Secondly, it is not always possible to uniquely detect two 
separate minima from the light curves of spotted stars. This is 
a general problem and common to all time series analysis meth- 
ods. Below some critical phase separation A^crft, the two minima 
merge into one common minimum. We calculated a theoretical 
limit of A0crit = 0.25 and then obtained A^cnt ~ 0.33 from sim- 
ulated data analysed with the CPS. The existence of this type 
of a phase limit is a concern when analysing photometry of ac- 
tive stars, since it means that two individual spots must have a 
considerable longitudinal separation to be detected separately. If 
one is not careful, this may lead to detections of spurious active 
longitudes. 

Thirdly, there is an intrinsic instability in the period estima- 
tion. This is a concern because period variations are often used to 
measure the differential rotation of active stars. For good quality 
data the effect is limited, but grows substantially for more noisy 
data. This period instability is caused by random effects due to 
noise, which is a general problem for all time series analysis 
methods, especially when the analysed datasets are short. 

We applied the CPS to 12 years of V band photometry of the 
young solar analogue HD 116956. The analysis revealed varia- 
tions in the mean magnitude M and the light curve amplitude A 
with an apparent period around 3 years. However, we found no 
conclusive evidence supporting the interpretation that these vari- 
ations are signs of an activity cycle, as the CPS analysis for the 
M and A estimates gave very high;^^^ values for a periodic fit. 

The star also displays two active longitudes that have re- 
mained stable over the whole 12 year observing period. The stan- 
dard Kuiper test of the primary and secondary minimum epochs 
gave an extreme critical level 2k = 8.7x10"^' for the rotation 
period P-^i = 7.8416 ± 0.001 1 d of these active longitudes. There 
have been only few transient excursions to other longitudes and 
the separation of the primary and secondary minima has stayed 
nearly unchanged. During the first observing season at 1999, 
the star underwent a flip-flop. In this event, the major activity 
switched from one active longitude to the other, while the old 
primary minimum disappeared completely. This flip-flop event 
happened very fast. 

We estimated the differential rotation of the star, assuming 
that it is the cause of the observed photometric period variations. 
These +3cr variations of the period gave the value Z - 11%, 
which is much larger than what could be caused only by spurious 
period variations. Assuming that the spot distribution covered 
the whole latitude range from equator to pole and that the solar 
law of differential rotation were valid also for HD 1 16956, these 
variations would correspond to a differential rotation coefficient 
\k\ = 0. 1 1 , or equivalently a differential rotation rate AH = 0.057 
rad d '. The observed solar value is A: = 0.20. If the latitudinal 
extent of the spot activity in HD 1 16956 were the same as in the 
Sun, surface differential rotation would be much stronger, i.e. 
\k\ = 0.44. 

The mean time scale of change of the light curve, i.e. the spot 
distribution of the star, was Tc = 44.1 d. This exceeds the length 
of the datasets, Ar^ax = 25.0 d. Hence the spot distribution re- 
mained unchanged during a typical dataset. There may be a link 
between the above Tc value and the convective turnover time 
Tc = 28.5 d, having the same order of magnitude. This could be 
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expected, if convective mixing causes the observed changes in 
the distribution of starspots. 
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Appendix A: Format of the results 

Through the CPS method we obtain a lot of information useful 
for studying stellar activity. These results are published electron- 
ically and it is therefore necessary to describe the notation for- 
mat in detail. The results of each individual dataset have been 
compressed to nine lines of an ASCII file. The parameter esti- 
mates are written in the format given in Table A. 1 . This table 
also specifies the units of the parameters. 

The first row in Table A. 1 gives the segment number (SEG), 
the epoch of the first observation in the segment (fo) and the 
dataset number (SET). The second row gives the epochs of the 
first (fi) and the last (f„) observation in the subset, the mean 
epoch of the dataset (t) and specifies whether the dataset is 
considered an independent dataset or not (IND). For indepen- 
dent datasets the value is 1ND=1, otherwise 1ND=0. The third 
row gives the number of observations (n), their mean (niy) and 
their standard deviation (sy). The fourth row gives the order of 
the model (K), the standard deviation of the residuals (cr^) and 
the time scale of change (Tc). As already mentioned in Sect. 
2.3, the value Tc - -2 indicates that the model describes well 
all the remaining datasets in the segment. The rest of the rows 
give values for the light curve parameters (M, P, A, /n,m,i, fmin,2) 
and their errors (ctm, o-p, cr^, o"(„,i„, , Cf„i„.2)- K is also specified 
whether these parameter estimates are considered reliable or not. 
Reliable estimates have R=0, while unreliable ones have R=l. 
All the heliocentric Julian dates are given in the truncated form 
HJD - 2 400 000. If for some reason no value has been obtained 
for some parameter, the "dummy" value - 1 is given. 

A typical example of an entry into the output file is: 



Table A.l. The format of the CPS output file for each individual datset. 





1 


51171.9835 




17 


5120Q, 


.9434 


51225.8738 


51214 


.2168 




22 


0.2658 


0, 


.0086 




2 


0.0055 


22, 


.2344 


0. 


,2654 


0.0012 







7. 


,1382 


0.1267 







0. 


,0274 


0.0045 







51207, 


.7204 


0.1689 







51203, 


.6602 


0.2811 




1 



This is the 17th dataset in the first segment (SEG=1, SET=17) 
in the analysis^ of the star HD 116956 presented in Sect. 4. It 
is a dataset of « = 22 observations obtained during a period of 
about Ar = f„ - fi = 25 d centered around HJD = 2451214.22. 
The dataset has been labelled as an independent one (IND=1). 
The modelling has been done with K - 2 and the values for all 
model parameters M, P, A, fmin.i and fn,in,2 have been obtained. 
Only the fmin,2 estimate is found to be unreliable (R = 1). 

Another example is the 11th dataset in the second segment 
(SEG=2, SET=11): 



[SEG] = integer 


[to] = HJD 


[SET] = integer 


[fi]=HJD 


[t„] = HJD 


[t] = HJD [IND] = integer 


[n] = integer 


[niy] = mag 


[Sy] = mag 


[K] = integer 


[o-J = mag 


[Tc] = d 


[M] = mag 


[cr„] = mag 


[R] = integer 


[P]=d 


[a-p] = d 


[R] = integer 


[A] = mag 


[cta] = mag 


[R] = integer 


[f„i„,i] = HJD 


['^v.„,]=d 


[R] = integer 


[f„,„,2] = HJD 


['^'^.nJ = d 


[R] = integer 


2 


51515.0410 


11 


51560.9421 


51585.8988 


51574.4281 


13 


0.2561 


0.0084 


1 


0.0063 


-1.0000 


0.2568 


0.0019 


1 


8.3355 


0.5855 


1 


0.0154 


0.0046 


1 


51567.3488 


0.7673 


1 


-1.0000 


-1.0000 


1 



This dataset is not independent (IND=0). It contains only n = 13 
observations and has been modelled with a simpler K - I model. 
For this model, there does not exist a secondary minimum and 
correspondingly fmin,2 has a dummy value -1. All light curve pa- 
rameters have been found unreliable (R = 1). Also, because the 
model parameters have been found unreliable, no value has been 
computed for the time scale of change Tc- 
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